library(yaml)
library(dplyr)
library(readr)
library(gridExtra)
library(ggplot2)
library(grid)
library(lspline)
library(lmtest)
library(car)
library(dfadjust)
CONFIG <- yaml.load_file('config_global.yaml')
build  <- CONFIG$build$descriptive
source(sprintf("%s/library.r", CONFIG$source$lib))

main <- function() {
  
  countries <- c("united_states", "switzerland", "france", "denmark", "canada", "new_zealand", "japan",
                 "australia", "britain", "norway", "sweden", "germany")
  formatted_names <- list("australia"         = "Australia",
                          "britain"           = "Britain",
                          "canada"            = "Canada",
                          "denmark"           = "Denmark",
                          "germany"    		    = "Germany",
                          "new_zealand"       = "New Zealand",
                          "norway"            = "Norway",
                          "sweden"            = "Sweden",
                          "switzerland"       = "Switzerland",
                          "united_states"     = "United States",
                          "united_states_alt" = "United States",
                          "switzerland_alt"   = "Switzerland",
                          "france"            = "France",
                          "japan"             = "Japan")
  
  all_data <- read.csv(sprintf("%s/data.csv", build))
  
  us_equality <- make_us_equality_table(all_data, countries)
  sd_table <- make_sd_benchmark_table(all_data, countries)
  
  # Plots
  BASE <- c()
  BIN  <- c()
  DATA <- c()
  LN   <- c()
  TOP2 <- c()
  SPLTpre <- c()
  SPLTpst <- c()
  IN      <- c()
  OUT     <- c()
  PVAL    <- c()
  for (c in c(countries, "united_states_alt", "switzerland_alt")){
    db     <- subset(all_data, country == c)
    
    if (c == "united_states"){
      m <- dfadjustSE(lm(partisan_affect_polarization_out ~ years, data = db))
      write_gslab_table(sprintf("%.03f", 2 * pt(-abs(m$coefficients["years", "Estimate"] / m$coefficients["years", "HC2 se"]), m$coefficients["years", "df"])), 
                        output = sprintf("%s/out_pvalue.txt", build),
                        header = "<tab:out_pvalue>")
    }
    
    # CSES Comparison
    if (c != "switzerland_alt" & c != "united_states_alt") make_cses_comparison_plot(db, all_data, c, formatted_names[[c]])
    
    # ANES - CSES comparison
    if (c == "united_states") make_anes_cses_plot(db, all_data)
    
    # East vs West Germany comparison
    if (c == "germany") make_east_west_germany_plot(db, all_data)
    
    # Base
    a <- plot_partisan_affect_polarization(db %>% select(years, partisan_affect_polarization,
                                                         question_wording, scaling),
                                           formatted_names[[c]], c, build, all_data)
    if (c != "switzerland_alt" & c != "united_states_alt") BASE <- cbind(BASE, a)
    
    # Bin
    a <- plot_partisan_affect_polarization(db %>% select(years, partisan_affect_polarization = partisan_affect_polarization_bin,
                                                         question_wording, scaling),
                                           formatted_names[[c]], sprintf("%s_bin", c), build, 
                                           all_data %>% select(years, country, partisan_affect_polarization = partisan_affect_polarization_bin,
                                                               question_wording))
    if (c != "switzerland_alt" & c != "united_states_alt") BIN <- cbind(BIN, a)
    
    # Data driven
    a <- plot_partisan_affect_polarization(db %>% select(years, partisan_affect_polarization = partisan_affect_polarization_data,
                                                         question_wording, scaling),
                                           formatted_names[[c]], sprintf("%s_data", c), build, 
                                           all_data %>% select(years, country, partisan_affect_polarization = partisan_affect_polarization_data,
                                                               question_wording))
    if (c != "switzerland_alt" & c != "united_states_alt") DATA <- cbind(DATA, a)
    
    # Leaner: Germany and Norway do not have a leaner question in any year
    #         Australia only has one year with the leaner question
    if (! c %in% c("australia", "germany", "norway")) {
      a <- plot_partisan_affect_polarization(db %>% 
                                               select(years, partisan_affect_polarization = partisan_affect_polarization_leaner,
                                                      question_wording, scaling, has_leaner)  %>%
                                               filter(has_leaner == T),
                                             formatted_names[[c]], sprintf("%s_leaner", c), build, 
                                             all_data %>% select(years, country, partisan_affect_polarization = partisan_affect_polarization_leaner,	question_wording))
      if (c != "switzerland_alt" & c != "united_states_alt") LN <- cbind(LN, a)
    } else {
      plot_no_leaners_question(c, formatted_names[[c]])
      if (c != "switzerland_alt" & c != "united_states_alt") LN <- cbind(LN, rbind(c("---   [---]"), c("(---, ---)")))
    }
    
    # Top 2
    db_longest_top_2 <- subset(db, as.character(top_2_parties) == getmode(as.character(top_2_parties)))
    a <- plot_partisan_affect_polarization(db_longest_top_2 %>% select(years, 
                                                                         partisan_affect_polarization = partisan_affect_polarization_top_two,
                                                                         question_wording, scaling),
                                             formatted_names[[c]], sprintf("%s_top_2", c), build,
                                             all_data %>% select(years, country,
                                                                 partisan_affect_polarization = partisan_affect_polarization_top_two,
                                                                 question_wording))
    if (c == "france"){
      #France only has 2 observations so SE is unidentified
      a[2] <- "(---, ---)"
    }
    if (c != "switzerland_alt" & c != "united_states_alt") TOP2 <- cbind(TOP2, a)

    # Split
    a <- plot_partisan_affect_polarization(db %>% select(years, partisan_affect_polarization,
                                                         question_wording, scaling),
                                           formatted_names[[c]], sprintf("%s_split", c), build, all_data, split = T)
    
    if (!grepl("_alt", c)) SPLTpre <- cbind(SPLTpre, a[1:2, , drop = F])
    if (!grepl("_alt", c)) SPLTpst <- cbind(SPLTpst, a[3:4, , drop = F])
    if (!grepl("_alt", c)) PVAL    <- cbind(PVAL, a[5, , drop = F])

    # In- vs Out-
    a <- plot_partisan_affect_polarization(db %>% select(years, partisan_affect_polarization_in, partisan_affect_polarization_out,
                                                         question_wording, scaling),
                                           formatted_names[[c]], sprintf("%s_in_out", c), build, all_data, in_out = T, YLAB = "Affect")
    if (c != "switzerland_alt" & c != "united_states_alt") IN  <- cbind(IN,  a[1:2, , drop = F])
    if (c != "switzerland_alt" & c != "united_states_alt") OUT <- cbind(OUT, a[3:4, , drop = F])
    
    # Supt 
    supt <- read.csv(sprintf("%s/sup_t_partisan_affect_polar_%s.csv", build, c), header = F)
    supt <- t(supt)
    colnames(supt) <- c("lower", "upper")
    temp <- cbind(db, supt)
    plot_partisan_affect_polarization(temp %>% select(years, partisan_affect_polarization, lower, upper,
                                                      question_wording, scaling),
                                      formatted_names[[c]], sprintf("%s_supt", c), build, all_data, supt = T)
  }
  
  # Make table
  TABLE <- rbind(BASE, us_equality$all$test_US_equality_p, us_equality$all$test_zero, TOP2, LN, DATA, BIN, 
                 SPLTpre, SPLTpst, us_equality$post$test_US_equality_p, us_equality$post$test_zero, PVAL,
                 sd_table[, "year"], sd_table[, "weighted_sd_of_pi"])
  write_gslab_table(TABLE, sprintf("%s/summary_table.txt", build), "<tab:summary_table>")
  
  # Make party share plots
  # and record the minimal share top 2 over survey years for each country
  min_share_top_2 <- matrix(NA, nrow = length(countries), ncol = 2)
  colnames(min_share_top_2) <- c("country", "min_share_top_2")
  for (i in c(1:length(countries))) {
    c     <- countries[i]
    top_2 <- read.csv(sprintf("build/descriptive/top_2_%s.csv", c))
    
    if (c %in% c("switzerland", "france", "japan")){
      top_2_cses <- read.csv("build/descriptive/top_2_cses.csv")
      if (c == "switzerland") top_2 <- rbind(top_2, top_2_cses[top_2_cses$year == "Switzerland_2003", ])
      if (c == "france") top_2      <- rbind(top_2, top_2_cses[grepl("France", top_2_cses$year), ])
      if (c == "japan") top_2       <- rbind(top_2, top_2_cses[grepl("Japan", top_2_cses$year), ])
      top_2$year <- as.numeric(gsub("Switzerland_|France_|Japan_", "", top_2$year))
    }
    
    plot_party_share(top_2, c, formatted_names[[c]], "share_top_2")
    plot_party_share(top_2, c, formatted_names[[c]], "share_has_party")
    min_share_top_2[i, "country"] <- c
    min_share_top_2[i, "min_share_top_2"] <- min(top_2$share_top_2)
  }
  write_gslab_table(min_share_top_2, sprintf("%s/min_share_top_2.txt", build), "<tab:min_share_top_2>")

}





make_us_equality_table <- function(all_data, countries){
  M <- c()
  S <- c()
  for (c in countries){
    m    <- dfadjustSE(lm(partisan_affect_polarization ~ years, subset(all_data, country == c)))
    s    <- dfadjustSE(lm(partisan_affect_polarization ~ lspline(years, c(2000)), subset(all_data, country == c)))
    
    if (c != 'united_states'){
      m_US <- dfadjustSE(lm(partisan_affect_polarization ~ years * as.numeric(country == 'united_states'), subset(all_data, country %in% c(c, 'united_states'))))
      s_US <- dfadjustSE(lm(partisan_affect_polarization ~ lspline(years, c(2000)) * as.numeric(country == 'united_states'), subset(all_data, country %in% c(c, 'united_states'))))
      
      M <- rbind(M, cbind(c, m$coefficients[2, 1], m$coefficients[2, 3], m$coefficients[2, 5], m_US$coefficients[4, 1], m_US$coefficients[4, 3], m_US$coefficients[4, 5]))
      S <- rbind(S, cbind(c, s$coefficients[3, 1], s$coefficients[3, 3], s$coefficients[3, 5], s_US$coefficients[6, 1], s_US$coefficients[6, 3], s_US$coefficients[6, 5]))
    } else {
      M <- rbind(M, cbind(c, m$coefficients[2, 1], m$coefficients[2, 3], m$coefficients[2, 5], 0, 1, 1))
      S <- rbind(S, cbind(c, s$coefficients[3, 1], s$coefficients[3, 3], s$coefficients[3, 5], 0, 1, 1))
    }
    
  }
  
  M           <- data.frame(M)
  colnames(M) <- c("country", "slope", "se", "dof", "slope_US", "se_US", "dof_US")
  M$slope     <- as.numeric(as.character(M$slope))
  M$se        <- as.numeric(as.character(M$se))
  M$slope_US  <- as.numeric(as.character(M$slope_US))
  M$se_US     <- as.numeric(as.character(M$se_US))
  M$dof       <- as.numeric(as.character(M$dof))
  M$dof_US    <- as.numeric(as.character(M$dof_US))
  
  S           <- data.frame(S)
  colnames(S) <- c("country", "slope", "se", "dof",  "slope_US", "se_US", "dof_US")
  S$slope     <- as.numeric(as.character(S$slope))
  S$se        <- as.numeric(as.character(S$se))
  S$slope_US  <- as.numeric(as.character(S$slope_US))
  S$se_US     <- as.numeric(as.character(S$se_US))
  S$dof       <- as.numeric(as.character(S$dof))
  S$dof_US    <- as.numeric(as.character(S$dof_US))
  
  M$test_US_equality_p <- 2 * pt(-abs(M$slope_US) / M$se_US, M$dof_US)
  S$test_US_equality_p <- 2 * pt(-abs(S$slope_US) / S$se_US, S$dof_US)
  
  M$test_zero          <- 2 * pt(-abs(M$slope) / M$se, M$dof)
  S$test_zero          <- 2 * pt(-abs(S$slope) / S$se, S$dof)
  
  return(list("all" = M, "post" = S))
}

make_sd_benchmark_table <- function(all_data, countries){
  b_table <- c()
  for (c in countries){
    temp <- subset(all_data, country == c)
    temp <- subset(temp, years == min(years))
    b_table <- rbind(b_table, 
                     cbind(c, min(temp$years), temp$affect_sample_sd))
  }
  colnames(b_table) <- c("country", "year", "weighted_sd_of_pi")
  return(b_table)
}

getmode <- function(v) {
  uniqv <- unique(v)
  return(uniqv[which.max(tabulate(match(v, uniqv)))])
}

plot_partisan_affect_polarization  <- function(db, name, extension, build, all_data, 
                                               YLAB = "Affective polarization", X_loc = 1970, supt = F, split = F, in_out = F){
  
  db$question_wording <- sprintf("%s %s", db$question_wording, db$scaling)
  
  if (in_out){
    in_  <- db[, c("years", "partisan_affect_polarization_in")]
    out_ <- db[, c("years", "partisan_affect_polarization_out")]
    out_$question_wording <- "Other"
    in_$question_wording  <- "Own"
    colnames(in_) <- colnames(out_) <- c("years", "partisan_affect_polarization", "question_wording")
    db <- rbind(in_, out_)
  }
  
  # https://stackoverflow.com/questions/7549694/adding-regression-line-equation-and-r2-on-graph
  # https://stackoverflow.com/questions/11618392/ggplot-text-printed-by-geom-text-is-not-clear
  lm_eqn <- function(df, table = F){
    m <- dfadjustSE(lm(partisan_affect_polarization ~ years, df))

    if (!table) return(sprintf("Slope: %.02f\n(%.02f, %.02f)", m$coefficients[2, 1], 
                               m$coefficients[2, 1] - 1.96 * m$coefficients[2, 4],
                               m$coefficients[2, 1] + 1.96 * m$coefficients[2, 4]))
    else {
      return(rbind(c(sprintf("%.02f %s [%.0f]", m$coefficients[2, 1], " ",
                             nrow(df))), 
                   c(sprintf("(%.02f, %.02f)",
                             m$coefficients[2, 1] - 1.96 * m$coefficients[2, 4],
                             m$coefficients[2, 1] + 1.96 * m$coefficients[2, 4]))))
    }
  }
  
  # Get Y Limits
  by_min <- ifelse(in_out, "partisan_affect_polarization_out", "partisan_affect_polarization")
  by_max <- ifelse(in_out, "partisan_affect_polarization_in",  "partisan_affect_polarization")
  c_min <- aggregate(all_data[, by_min],
                     by = list(all_data$country), FUN = min, na.rm = T)
  c_max <- aggregate(all_data[, by_max],
                     by = list(all_data$country), FUN = max, na.rm = T)
  colnames(c_min) <- c("country", "min")
  colnames(c_max) <- c("country", "max")
  c_       <- merge(c_min, c_max)
  diff_max <- max(c_$max - c_$min)
  Ymin <- min(db$partisan_affect_polarization)
  Ymax <- max(db$partisan_affect_polarization)
  YLIM_max <- Ymin + (Ymax - Ymin) / 2 + diff_max / 2 + 14 + (supt * 4) + (in_out * 18)
  YLIM_min <- Ymin + (Ymax - Ymin) / 2 - diff_max / 2 - (supt * 4) - (in_out * 15)  - 1
  
  # Replace US
  if (name == "United States"){
    name <- "US"
  }
  
  # Plot
  y_by <- ifelse(in_out, 10, 5)
  plt  <-  ggplot(db, aes(x = years, y = partisan_affect_polarization)) +
    labs(title = name, x = "", y = ifelse(grepl("united_states|switzerland_alt", extension), YLAB, "")) +
    xlim(min(all_data$years), max(all_data$years)) +
    scale_y_continuous(limits = c(YLIM_min, YLIM_max),
                       breaks = seq(ceiling(YLIM_min / y_by) * y_by, floor(YLIM_max / y_by) * y_by, y_by))
  
  
  if (supt){
    plt <- plt + geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) 
  }
  
  # Standardize shape assignment of question wording across plots
  db$question_wording <- factor(db$question_wording)
  old.lvl<-levels(db$question_wording)
  db$question_wording <- factor(db$question_wording,
                                levels=c("favourable (0-100)", "like (0-10)",
                                         "like (0-100)",       "Own",
                                         sort(old.lvl[old.lvl   != "like (0-10)"
                                                      & old.lvl != "favourable (0-100)"
                                                      & old.lvl != "like (0-100)"
                                                      & old.lvl != "Own"], 
                                              decreasing=T)))
  
  if (split){
    pre_2000  <- subset(db, years <= 2000)
    post_2000 <- subset(db, years > 2000)
    
    out  <- lm(partisan_affect_polarization ~ lspline(years, c(2000)), data = db)
    lab1 <- sprintf("Slope: %.02f\n(%.02f, %.02f)", dfadjustSE(out)$coefficients[2, 1], 
                    dfadjustSE(out)$coefficients[2, 1] - 1.96 * dfadjustSE(out)$coefficients[2, 4], 
                    dfadjustSE(out)$coefficients[2, 1] + 1.96 * dfadjustSE(out)$coefficients[2, 4])
    lab2 <- sprintf("Slope: %.02f\n(%.02f, %.02f)", dfadjustSE(out)$coefficients[3, 1], 
                    dfadjustSE(out)$coefficients[3, 1] - 1.96 * dfadjustSE(out)$coefficients[3, 4],
                    dfadjustSE(out)$coefficients[3, 1] + 1.96 * dfadjustSE(out)$coefficients[3, 4])
    
    plt <- plt + geom_point(aes(shape = db$question_wording), show.legend = F, size = 2.2) +
      geom_segment(aes(x = min(db$years), y = predict(out, newdata = data.frame(years = c(min(db$years)))), 
                       xend = 2000, yend = predict(out, newdata = data.frame(years = c(2000)))), color = "black", show.legend = F) +
      geom_segment(aes(x = 2000), y = predict(out, newdata = data.frame(years = c(2000))), 
                       xend = max(db$years), yend = predict(out, newdata = data.frame(years = c(max(db$years)))), color = "red4", show.legend = F) + 
      annotate("text", x = X_loc,  y = (YLIM_max - YLIM_min) * .95 + YLIM_min,
               label = lab1,  colour = "black", hjust = 0, size = 4.5) +
      annotate("text", x = 2013,   y = (YLIM_max - YLIM_min) * .95 + YLIM_min, 
               label = lab2, colour = "red4",  hjust = 1, size = 4.5)
    
  } else if (in_out){
    plt <- plt + geom_point(aes(color = db$question_wording, shape = db$question_wording), show.legend = T, size = 2.2) +
      scale_colour_manual(name = "Party",
                          values = c("black", "red4"), 
                          labels = c("Own", "Other")) +
      scale_shape_manual(name = "Party", 
                         values = c(19, 17),
                         labels = c("Own", "Other")) +
      geom_smooth(data = in_, color = "black",
                  method = "lm", formula = y~x, se = F, show.legend = F) +
      geom_smooth(data = out_, color = "red4",
                  method = "lm", formula = y~x, se = F, show.legend = F) +
      annotate("text", x = X_loc,  y = (YLIM_max - YLIM_min) * .95 + YLIM_min,
               label = lm_eqn(in_),  colour = "black", hjust = 0, size = 4.5) +
      annotate("text", x = X_loc,   y = (YLIM_max - YLIM_min) * .05 + YLIM_min,
               label = lm_eqn(out_), colour = "red4",  hjust = 0, size = 4.5) 
  } else {
    plt <- plt + geom_point(aes(shape = db$question_wording), size = 2.2) +
      geom_smooth(color = "red4", method = "lm", formula = y~x, se = F, show.legend = F) + 
      annotate("text", x = X_loc, y = (YLIM_max - YLIM_min) * .95 + YLIM_min,
               label = lm_eqn(db), colour = "red4", hjust = 0, size = 4.5) +
      scale_shape_discrete(name = "Original response scaling")
  }
  
  # Edit Theme
  plt <- plt + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                                  panel.grid.minor = element_blank(),
                                  axis.line = element_line(colour = "black"),
                                  legend.position = c(.75, .83), 
                                  legend.text  = element_text(size = 11), 
                                  legend.title = element_text(size = 11), 
                                  plot.title   = element_text(hjust = 0.5, size = 20),
                                  axis.title.x = element_text(size = 16),
                                  axis.title.y = element_text(size = 16),
                                  axis.text.x  = element_text(size = 16),
                                  axis.text.y  = element_text(size = 16))
  
  ggsave(sprintf("%s/%s_partisan_affect_polarization_plot.pdf", build,
                 tolower(extension)), height = 4.2, width = 5.5)
  
  if (split){
    out_pre_post_test  <- dfadjustSE(lm(partisan_affect_polarization ~ lspline(years, c(2000), marginal = T), data = db))

    table_out <- rbind(c(sprintf("%.02f %s [%.0f]", dfadjustSE(out)$coefficients[2, 1], " ",
                                 sum(db$years <= 2000))), 
                       c(sprintf("(%.02f, %.02f)", 
                                 dfadjustSE(out)$coefficients[2, 1] - 1.96 * dfadjustSE(out)$coefficients[2, 4], 
                                 dfadjustSE(out)$coefficients[2, 1] + 1.96 * dfadjustSE(out)$coefficients[2, 4])))
    table_out <- rbind(table_out,
                       c(sprintf("%.02f %s [%.0f]", dfadjustSE(out)$coefficients[3, 1], " ",
                                 sum(db$years > 2000))), 
                       c(sprintf("(%.02f, %.02f)", 
                                 dfadjustSE(out)$coefficients[3, 1] - 1.96 * dfadjustSE(out)$coefficients[3, 4], 
                                 dfadjustSE(out)$coefficients[3, 1] + 1.96 * dfadjustSE(out)$coefficients[3, 4])))
    
    table_out <- rbind(table_out, sprintf("%.03f", 2 * pt(-abs(out_pre_post_test$coefficients[3, 1]) / out_pre_post_test$coefficients[3, 3], out_pre_post_test$coefficients[3, 5])))
  } else if (in_out){
    table_out <- lm_eqn(in_, table = T)
    table_out <- rbind(table_out, lm_eqn(out_, table = T))
  } else {
    table_out <- lm_eqn(db, table = T)
  }
  return(table_out)
}

plot_party_share <- function(db, c, name, series){
  # https://stackoverflow.com/questions/10349206/add-legend-to-ggplot2-line-plot
  db$group <- sprintf("%s:%s", db$top_party_1, db$top_party_2)
  
  db           <- db[, c("year", "group", series)]
  colnames(db) <- c("year", "group", "share")
  
  plt <- ggplot(db, show.legend = F) + xlim(1960, 2020) +
    scale_y_continuous(limits = c(0, 1)) + geom_line(aes(x = year, y = share), show.legend = F) +
    scale_colour_manual(values=c("black","red4"))
  
  plt <- plt + labs(title = name, x = "", y = ifelse(name %in% c("US", "United States"), "Proportion", ""), color = "")
  if (series == "share_top_2"){
    plt <- plt + geom_point(aes(x = year, y = share, shape = factor(group)), size = 2.2) + scale_shape_discrete(guide = FALSE)
  } else {
    plt <- plt + geom_point(aes(x = year, y = share), size = 2.2) 
  }
  
  # Edit Theme
  plt <- plt + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                                  panel.grid.minor = element_blank(),
                                  axis.line = element_line(colour = "black"),
                                  plot.title   = element_text(hjust = 0.5, size = 18),
                                  axis.title.x = element_text(size = 14),
                                  axis.title.y = element_text(size = 14),
                                  axis.text.x  = element_text(size = 12),
                                  axis.text.y  = element_text(size = 12))
  
  ggsave(sprintf("%s/%s_%s_plot.pdf", build, c, series), height = 3.2, width = 5.5)
}

plot_no_leaners_question <- function(c, name) {
  TEXT <- ifelse(c == "australia", "Only one survey\n with leaners question", "No leaners question")
  pdf(sprintf('%s/%s_leaner_partisan_affect_polarization_plot.pdf', build, c), height = 4.2, width = 5.5)
  par(mar = c(0,0,0,0))
  plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
  text(x = 0.5, y = 0.5, TEXT, cex = 1.6, col = "black")
  text(x = 0.5, y = 0.99, sprintf("%s", name), cex = 1.6, col = "black")
  dev.off()
}

make_anes_cses_plot <- function(db, all_data){

  # Load Data and Clean
  boot <- read.csv("build/descriptive/anes_cses_bootstrap.txt",  header = F, skip = 1, stringsAsFactors = F)
  data <- read.csv("build/descriptive/anes_cses_comparison.txt", header = F)
  
  ANES           <- data[, c(2, 1)]
  CSES           <- data[, c(3, 1)]
  colnames(ANES) <- colnames(CSES) <- c("Affect", "Year")
  ANES$Type      <- "ANES: 0-100 (feel)"
  CSES$Type      <- "CSES: 0-10 (like)"
  
  data      <- rbind(ANES, CSES)
  data$Type <- factor(data$Type)
  
  
  # Plot
  two_line_plot(db, limit_data = all_data, plotting_data = data, Black = ANES, Red = CSES,
                name = "anes_cses_comparison",
                black_label = "ANES: 0-100 (feel)", red_label = "CSES: 0-10 (like)", 
                black_text_loc = c(1997, 32), red_text_loc = c(1996, 45), difference = T, country = "US", boot = boot)
  

}

make_cses_comparison_plot <- function(db, all_data, country, country_name){
  
  orig_countries <- c("australia", "britain",    "canada", "germany", "new_zealand", "norway", "sweden", "switzerland", "united_states", "denmark", "france", "japan")
  cses_countries <- c("Australia_", "Great Britain_", "Canada_", "Germany_", "New Zealand_", "Norway_", "Sweden_", "Switzerland_", "United States of America_", "Denmark_", "France_", "Japan_")
  cses_country   <- cses_countries[which(orig_countries == country)]
  
  # Load Data and Clean
  cses <- read.csv(sprintf("%s/cses_only.csv", build), stringsAsFactors = F)

  cses <- cses %>% filter(grepl(cses_country, years))
  cses$years <- as.numeric(gsub(cses_country, "", cses$years))
  first_cses_year <- min(cses$years)
  if (sum(db$years %in% c(first_cses_year, first_cses_year - 1)) == 1){
    cses$partisan_affect_polarization <- cses$partisan_affect_polarization + (db$partisan_affect_polarization[db$years %in% c(first_cses_year, first_cses_year - 1)] - cses$partisan_affect_polarization[cses$years == first_cses_year])
  } else {
    cses$partisan_affect_polarization <- cses$partisan_affect_polarization + (db$partisan_affect_polarization[db$years %in% c(first_cses_year)] - cses$partisan_affect_polarization[cses$years == first_cses_year])
  }
  
  cses$Affect <- cses$partisan_affect_polarization
  db$Affect   <- db$partisan_affect_polarization
  cses$Year   <- cses$years
  db$Year     <- db$years
  cses$Type   <- "CSES"
  db$Type     <- "Our main series"
  
  cses$Year <- cses$Year - 1/3
  
  data      <- rbind(db, cses)
  data$Type <- factor(data$Type)  
  
  if (sum(db$years %in% c(first_cses_year, first_cses_year - 1)) == 1){
    db <- db %>% filter(Year >= (first_cses_year - 1))
  } else {
    db <- db %>% filter(Year >= first_cses_year)
  }
   
  # Plot
  two_line_plot(db, limit_data = all_data, plotting_data = data, Black = cses, Red = db,
                name = sprintf("%s_cses_comparison", country),
                black_label = "CSES", red_label = "Our Main Series", 
                black_text_loc = c(1977, 35), red_text_loc = c(1977, 60), difference = F, country = country_name, lines = F, scale = 1, 
                YLAB = ifelse(country == "united_states", "Affective polarization", ""), ylim_drop = T)
  
}

two_line_plot <- function(db, limit_data, plotting_data, Black, Red, name, black_label, red_label, black_text_loc, red_text_loc, difference, country, 
                          boot = NULL, lines = T, scale = 1.5, YLAB = "Affective polarization", ylim_drop = F){
  
  # Temp Helper Function
  lm_eqn <- function(df, boot){
    m <- dfadjustSE(lm(Affect ~ Year, df))
    if (is.null(boot)){
      return(sprintf("Slope: %.02f\n(%.02f,%.02f)", m$coefficients[2, 1], 
                     m$coefficients[2, 1] - 1.96 * m$coefficients[2, 4],
                     m$coefficients[2, 1] + 1.96 * m$coefficients[2, 4]))
    } else {
      return(sprintf("Slope: %.02f", m$coefficients[2, 1]))
    }
  }
  
  # Get Y Limits
  y_by <- 5
  c_min <- aggregate(limit_data[, "partisan_affect_polarization"],
                     by = list(limit_data$country), FUN = min, na.rm = T)
  c_max <- aggregate(limit_data[, "partisan_affect_polarization"],
                     by = list(limit_data$country), FUN = max, na.rm = T)
  colnames(c_min) <- c("country", "min")
  colnames(c_max) <- c("country", "max")
  c_       <- merge(c_min, c_max)
  diff_max <- max(c_$max - c_$min)
  Ymin <- min(db$partisan_affect_polarization)
  Ymax <- max(db$partisan_affect_polarization)
  YLIM_max <- Ymin + (Ymax - Ymin) / 2 + diff_max / 2 + 10 + 4
  YLIM_min <- Ymin + (Ymax - Ymin) / 2 - diff_max / 2 - 4 - (ylim_drop * 5)
  
  # Plot (https://stackoverflow.com/questions/12410908/combine-legends-for-color-and-shape-into-a-single-legend)
  plt <- ggplot(plotting_data, aes(x = Year, y = Affect))  +
    xlim(min(limit_data$years), max(limit_data$years)) +
    scale_y_continuous(limits = c(YLIM_min, YLIM_max),
                       breaks = seq(ceiling(YLIM_min / y_by) * y_by, floor(YLIM_max / y_by) * y_by, y_by)) + 
    labs(title = country, x = "", y = YLAB)

  plt <- plt + geom_point(aes(color = Type, shape = Type), show.legend = T, size = 2.2) +
    scale_colour_manual(name = "Type",
                        values = c("black", "red4"), 
                        labels = c(black_label, red_label)) +
    scale_shape_manual(name = "Type", 
                       values = c(19, 17),
                       labels = c(black_label, red_label))
    
  if (lines){
    plt <- plt + geom_smooth(data = Black, color = "black",
                                method = "lm", formula = y~x, se = F, show.legend = F) +
        geom_smooth(data = Red, color = "red4",
                    method = "lm", formula = y~x, se = F, show.legend = F) +
        annotate("text", x = black_text_loc[1], y = black_text_loc[2],
                 label = lm_eqn(Black, boot),  colour = "black", hjust = 0, size = 4.5) +
        annotate("text", x = red_text_loc[1], y = red_text_loc[2],
                 label = lm_eqn(Red, boot), colour = "red4",  hjust = 0, size = 4.5)
  }
   
  if (difference){
    plt <- plt + annotate("text", x = 1980, y = 37, 
                          label = sprintf("Difference: %.02f\nSE: %.02f", as.numeric(boot[1, 1]), as.numeric(boot[2, 1])),
                          hjust = 0, size = 4.5, colour = "grey30")
  }

  # Edit Theme
  plt <- plt + theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                                  panel.grid.minor = element_blank(),
                                  axis.line = element_line(colour = "black"),
                                  legend.position = c(.8, .85), 
                                  legend.text  = element_text(size = 11), 
                                  legend.title = element_text(size = 11), 
                                  plot.title   = element_text(hjust = 0.5, size = 20),
                                  axis.title.x = element_text(size = 16),
                                  axis.title.y = element_text(size = 16),
                                  axis.text.x  = element_text(size = 16),
                                  axis.text.y  = element_text(size = 16))
  
  ggsave(sprintf("%s/%s.pdf", build, name), height = 4.2 * scale, width = 5.5 * scale)
}

make_east_west_germany_plot <- function(db, all_data){
  # Load Data and Clean
  data <- all_data %>% 
    filter(country == "germany" | country == "germany_alt") %>%
    filter(!years %in% c(1996, 1997, 1998)) %>%
    select(partisan_affect_polarization, years, country) %>%
    mutate(country = if_else(country == "germany_alt", "East: -5 to 5 (think highly)", "West: -5 to 5 (think highly)"))
  
  colnames(data) <- c("Affect", "Year", "Region")
  East <- data %>% filter(Region == "East: -5 to 5 (think highly)")
  West <- data %>% filter(Region == "West: -5 to 5 (think highly)")

  East$Type      <- "East: -5 to 5 (think highly)"
  West$Type      <- "West: -5 to 5 (think highly)"
  
  data      <- rbind(East, West)
  data$Type <- factor(data$Type)  
  
  # Plot
  two_line_plot(db, limit_data = all_data, plotting_data = data, Black = East, Red = West, name = "east_west_germany_comparison",
                black_label = "East: -5 to 5 (think highly)", red_label = "West: -5 to 5 (think highly)", 
                black_text_loc = c(2010, 40), red_text_loc = c(1980, 30), country = "Germany", difference = F)
}

main()